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Oh 

O An imperative condition for the functioning of a power-grid network is that its power gener- 

ators remain synchronized. Disturbances can prompt desynchronization, which is a process 

Q that has been involved in large power outages. Here we derive a condition under which the 

^ desired synchronous state of a power grid is stable, and use this condition to identify tunable 

4^ parameters of the generators that are determinants of spontaneous synchronization. Our 
Q-i 

I— I analysis gives rise to an approach to specify parameter assignments that can enhance syn- 
chronization of any given network, which we demonstrate for a selection of both test systems 

t:;^- and real power grids. Because our results concern spontaneous synchronization, they are 

^ relevant both for reducing dependence on conventional control devices, thus offering an ad- 

^ ditional layer of protection given that most power outages involve equipment or operational 

^ errors, and for contributing to the development of "smart grids" that can recover from fail- 

^ ures in real time. 

> 
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The current resounding interest in network synchronizatioiJiEl has been stimulated by the 
prospect that theoretical studies will help explain the behavior of real complex network^. Re- 
cent advances include the modeling of chimera spatiotemporal pattern^, the discovery of low- 
dimensional dynamics in heterogeneous population^, and the characterization of network syn- 
chronization landscape^. However, despite the significant insights provided by these and other 
studiePI^, the connection between network theory and the synchronization of real systems re- 
mains under-explored. This is partly due to the scarcity of dynamical data, partly due to the ide- 
alized nature of the theoretical constructs (for network-unrelated exceptions, see refs. 16 and[T7|). 



In light of these considerations, we explore power-grid systems as genuine complex networks of 
broad significance that are amenable to theoretical modeling and whose dynamics can be simu- 
lated reliably. As the largest man-made machines in existence, modem power-grid networks often 
consist of thousands of power substations and generators linked across thousands of kilometers. 
This complexity is reflected in a rich variety of collective behaviors and instabilities these systems 
reportedly exhibitpl. 

Here, we focus on the synchronization of power generators — a fascinating phenomenon that 
can occur spontaneously, that keeps all connected generators in pace, and whose failure constitutes 
an important source of instabilities in power-grid systems. In a network of n alternating current 
generators, a synchronous state is characterized by 

Si = S2 = ■ ■ ■ = 6n, (1) 

where 5i = 5i{t) represents the rotational phase of the ith generator and the dot represents the 
time derivative. This synchronization frequency is assumed to be close (albeit not necessarily 
equal) to a reference frequency. That coupled power generators can synchronize spontaneously 
is well known, as popularized in ref. \19\ and this has generated recent interest in the physics 
communit5E2ll2ll. However, the governing factors and hence the extent to which this phenomenon 
may occur in real power grids remain essentially unaddressed. Among the other studies that have 
provided substantive insights into power grid synchronization we mention the analysis of nonlin- 
ear modes associated with instabilities (refs. 22-|24 and references therein) and of an equivalence 



between power systems and networks of nonuniform Kuramoto oscillatorp^^. The former stud- 
ies provide an informative decomposition of nonlinear oscillations that arises when the network 
loses synchrony. The latter studies establish a sufficient condition for a power grid to converge to a 
synchronous state, which is obtained via singular perturbation analysis from a similar condition for 
nonuniform Kuramoto oscillators^^j^. This condition is necessary and sufficient for non-identical 
dynamical units under the assumption that the coupling strengths are unifornP^. In contrast, the 
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stability condition we derive below is necessary and sufficient for heterogeneously coupled power 
grids, with the assumption that a certain function of generator parameters is homogeneous. This 
condition helps us address the question of how to strengthen synchronization. 

Power grids deliver a growing share of the energy consumed in the world and will soon un- 
dergo substantial changes owing to the increased harnessing of intermittent energy sources, the 
commercialization of plug-in electric automobiles, and the development of real-time pricing and 
two-way energy exchange technologies. These advances will further increase the economical and 
societal importance of power grids, but they will also lead to new disturbances associated with 
fluctuations in production and demand, which may trigger desynchronization of power generators. 
Although power grids can, and often do, rely on active control devices to maintain synchronism, 
the question of whether proper design would allow the same to be achieved while reducing de- 
pendence on existing controllers is extremely relevant. In the U.S., for example, data on reported 
power outages show that over 3/4 of all large events involve equipment misoperation or human 
errors among other factord^. This illustrates why stability drawn from the network itself would be 
desirable. 



The dynamics 

We represent the power grid as a network of power generators and substations (nodes) connected 
by power transmission lines (links). The nodes may thus consume, produce, and distribute power, 
while the links transport power and may include passive elements with resistance, capacitance, and 
inductance. The state of the system is determined using power flow calculations given the power 
demand and other properties of the system (see Methods), which is a procedure we implement 
for the systems in Table 1. In possession of all variables that describe the steady-state alternating 
current flow, we seek to identify the conditions under which the generators can remain stably 
synchronized. 

The starting point of our analysis is the equation of motion 

2Hi (f5i 



-Pmi -Peii (2) 



the so-called swing equation (ref. 29 and Methods), which, along with the power flow equations, 
describes the dynamics of generator i. The parameter Hi is the inertia constant of the generator, cor 
is the reference frequency of the system, P^j is the mechanical power provided by the generator, 
and Pei is the power demanded of the generator by the network (including the power lost to damp- 
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ing). In equilibrium, P^i = Pei and the frequency coi = 5i remains equal to a common constant for 
all i, which is a sufficient condition for the frequency synchronization expressed in equation ([T]). 

When the system changes because of a fault or a fluctuation in power demand, so does Pei- 
According to equation ([2]), immediately after the change the difference between the power de- 
manded and the power supplied by the generator is compensated by either increasing or decreas- 
ing the angular momentum from the generator's rotor. The change in angular momentum does not 
remove the difference between Pei and Pmi', it only compensates for it until the spontaneous or 
controlled response of Pgj and P^j can adjust the r.h.s. of equation ([2]). It is through this process 
that the generators can re-achieve the synchrony they may have lost in entering the transient state. 
But what are the properties of the network connecting the generators? 

To address this question, we will first represent the power consumed at nodes of the network 
as equivalent impedance^. This involves the assumption that the network structure is fixed and 
the power demand is constant, which is appropriate because we are interested in the stability of 
the system on short time scales. For a given component, the corresponding admittance Y (the 
reciprocal of the impedance) is related to the voltage difference V and the complex conjugate of 
the power P through F = P/|l^p. Thus, the admittances are usually complex. 

The network structure 

The admittance matrix Yo= (^oij) is defined analogously to the Laplacian matrix in oscillator 
networks: loii is the negative of the admittance between nodes i and j ^ i, while Yqh is the sum 
of all admittances connected to node i (the row sum is nonzero, however, due to finite impedances 
to the ground). If V is the vector of node voltages and I is the vector of currents injected into the 
system at the nodes, then, according to the Kirchhoff's law, I = YqV. Because the non-generator 
nodes are now regarded as constant impedances, all nodes have zero injection currents except for 
the generator nodes; this can be used to represent the effective interactions among the generators 
through a reduced network. 

Indeed, rearranging Yq such that the first n indices correspond to the generator nodes and the 
remaining r indices correspond to non-generator nodes, one has 
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The system is then converted to I„ = YV„ by eliminating V^, a procedure known as Kron 
reductioil^, where the resulting effective admittance matrix is 

Y = Yonxn YonxrYo^xrYorxn- (4) 

The existence of the matrix Yo~x,,. follows from the assumed uniqueness of the voltage vector V 
(with respect to a reference voltage). In a system with n generators, the symmetric n x n matrix 
Y defines a network in which generator nodes i and j ^ i are connected by a link of strength 
determined by the effective admittance —Yij. The effective admittances are dominantly imaginary 
(the real part is smaller by up to one order of magnitude in the power grids we consider, as shown 
in Table 2). This indicates, as it will become clear below, that capacitances and inductances play 
an important role in the synchronization of generators. 

The effective network connecting the generator nodes, represented by Y, is significantly 
different from the physical network of transmission lines represented by Yq, as shown in Fig. 1 for 
the power grid of Northern Italy. If the physical network is connected, as observed for this and 
any other system that we can treat as a single power grid, then the effective network will be fully 
connected!^, meaning that any two generators are directly coupled through an effective admittance. 
This property of the effective network holds true under the assumption that the network is in a 
steady state, and is in contrast to the structure of the corresponding physical network, which is 
typically only sparsely connected. The connection strengths are not uniform in both the physical 
and the effective network (Fig. 1). However, the distribution of weighted degrees (the sum of the 
absolute values of all admittances connected to a node) is generally more homogeneous for Y than 
for Yo, as illustrated in Table 1 for several power networks. This degree homogeneity is a property 
of interest since it is known that it facilitates synchronization in networks of diffusively coupled 
oscillatord^^M]. Surprisingly, as we show next, in power-grid networks this property can be less 
determinant for synchronization than is the relation of certain generator parameters to the network 
structure. 

The stability condition 

To establish an interpretable relationship between synchronization robustness and power-grid pa- 
rameters, as well as to provide a condition that is necessary for any stricter form of stability, 
we focus on the short-term stability of the synchronous states under small perturbations. We 
thus linearize equation ([2]) around an equilibrium (synchronous) state, associated with electrical 
power Pg* and mechanical power P^^ and represented by S* and to*. Assuming that Si = S* + S'^, 
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Pei = P*i + Pgj, and Pmi = P^i + P^i, the linearized equation reads 



2H, d?5' 



dP dP 



duji 



(5) 



where we neglect the dependence of the mechanical power on changes in the phase Si, and denote 
u'^ = 5^. For the first term on the right-hand side, we assume the droop equation dPmi/dui = 
— l/{ujjiRi), where i?j > is a regulation parameter. For the second term, we assume there is a 
constant damping coefficient Di > such that dPei/dcOi = Di/ur. The third term can be obtained 
from P^ = DiJJuR + ^^=1 EiE^^Bij cos 5*^ - sin 5*^)5'^^, where 5*^ = 6* - 6*, =6[- 5'j, 
Gij and Bij are the real and imaginary components of Yij, and is the internal- voltage magnitude 
of the ith generatoi^^. We denote the n-dimensional vectors of 5^ and 6^ by Xi and X2, respectively. 



This yields a coupled set of 2n first-order equations, 

Xi = X2, 

X2 = -PX1-BX2, 

where P = (Pij) is the zero row sum matrix given by 



(6) 
(7) 



ujREiE. 



P 



(8) 



and B is the diagonal matrix of elements /3i = {Di + 1/ Ri) /2Hi. If we now assume this parameter 
/3i to be the same for all generator nodes, equations ([6]) and (|7]) can be diagonalized using the 
substitution Zi = Q^^Xi and Z2 = Q^^X2, where Q is a matrix of eigenvectors of P so that J 
= Q~^PQ is the diagonal matrix of the corresponding eigenvalues. This leads to Zi = Z2 and 
Z2 = — JZi — /3Z2, where (3 is the common value of the parameter Naturally, (3i is not the 
same for all generators in the power grids we consider, as evidenced by the normalized standard 
deviation in Table 1 (column 7), but this auxiliary assumption will allow us to derive results to 
enhance synchronization stability of power grids with nonidentical (3i. 

Equations ([6]) and (|7]) are then reduced to n decoupled 2-dimensional systems of the form 
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where aj is the jth eigenvalue of the coupling matrix P. The matrix P always has a null eigenvalue 
associated with a uniform shift of all phases, determined by Xi oc (1, . . . , 1)^, which corresponds 
to two perturbation eigenmodes (one for X2 = and one for X2 = — /3Xi) that do not affect 
the synchronization condition in equation ([T]). We denote this eigenvalue by ai and exclude the 
corresponding eigenmodes from subsequent analysis. 

The stability of the synchronous state is governed by the Lyapunov exponents of system (|9]), 

/3 , 1 



Xj^{a„P) = -| ± -v//3'-4«,-. (10) 

The synchronous state is stable iff the real component of every Lyapunov exponent is negative for 
all j > 2. That is, the stability condition translates to 

max ReXj^ < 0, for j = 2, . . . , n. (1 1) 



Because equation (|9]) is formally identical for all j, it is useful to drop the index j in equation ( fTO| ) 



and define the function A/3(a;) = max{-|-} ReX±{a, (3) for fixed (3 and tunable parameter a. This 
function can be interpreted as a master stability function: the stability of the system is determined 
by whether all the eigenvalues 02, • • • ,«n fall within the negative region of A^. This stability 



condition is of the form of those in refs. 34 and 35 for coupled oscillators, as explicitly shown in 
Methods. 

The matrix P is generally asymmetric, meaning that the influence of a generator on another 
generator is not necessarily identical to the influence of the second on the first. But we argue that 
the eigenvalues of P are all real over a range of conditions. We first factor the matrix as P = H~^P', 
where H is the diagonal matrix of inertia constants Hi, to then note (inspired by a similar ar- 
gument in ref. |33j) that the eigenvalues of P equal the eigenvalues of the partially symmetrized 
matrix P" = H~^/^P'H^^/^. The matrix P" is a zero row sum matrix of off-diagonal elements 
^rE^Ej sin 5*- — Bij cos 6* A. We now observe that, due to the small real-to-imaginary com- 
ponent ratio of the admittances and the predominance of small phase differences 5*j (Table 2), the 
antisymmetric part of P" is usually much smaller than the symmetric one; the 2-norm is over one 
order of magnitude smaller for the systems considered here, as shown in Table 3. If the eigenvalues 
of the symmetric part of P", which are necessarily real, are nondegenerate (as generally observed 
in practice), then we can show that the eigenvalues remain strictly real in the presence of a small 
antisymmetric component. Therefore, while our stability condition can be used in general, in what 
follows we will assume that P has real eigenvalue spectra, a condition observed to be satisfied in 
all of our simulations. 
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For real a, the stability region defined by Kp{a) < corresponds to a > 0. Indeed, the 
function is zero for a = and becomes positive and increasingly larger for smaller a < 0; 
in the interval < a < /3^/4, the function is negative and decreases as a increases. For 
a > 13'^/ 4:, the function A^ is negative and remains constant since the radicand in equation ([10]) is 
negative, rendering to be complex. This behavior is summarized in Fig. 2a for different values 
of /3, which further illustrates that varying /3 changes the behavior of the Lyapunov exponents as 
a function of a and, in particular, A^(a > /3^/4) = —(3/2, but the stability condition that aj > 
for all j > 2 does not change. Moreover, for a given (3 > 0, the master stability condition does not 
depend on whether this factor (3 is due to the presence of the regulation parameter Ri, the presence 
of damping Di, or a combination of both. However, all other parameters left constant, both (3 and 
the eigenvalues aj decrease as the inertia parameters Hi are uniformly increased. This has the 
effect of reducing |_ReAj^ | and hence the strength of the linear stability or instability. 



Enhancement of synchronization stability 

We now apply this formalism to enhance the stability of the synchronous states by tuning the 
parameters of the generators. For (3i = (3, the synchronization stability is solely determined by 
the master stability function Afs{a) at a = a2, the smallest nonzero (real) eigenvalue of P. This 
follows from the non- increasing dependence of A/3 (a) on a. For 02 > 0, which corresponds to 
stable synchronization, the quantity A^(a2) as a function of /3 attains its minimum at 

/9 = /Sopt = (12) 

(see Fig. 2b). Given a network structure and a steady state, and hence a fixed value of 0:2, this point 
of maximum stability can be achieved by adjusting one of the generator parameters. 

Specifically, one can ensure that /3j = /3opt for all generators by adjusting the droop parameter 

Ri to 

Ri 

or the damping coefficient Di to 

while keeping the other parameters unchanged. The tuning of the former, Ri, is most suitable for 
off-line optimization of stability, since the time scales associated with this parametei^^^ are usually 
larger than those associated with typical instabilities. The latter, D^, can be adjusted dynamically at 



4Hi^ - Di 



z = 1, . . . 



(13) 



4:Hiy/a^ - —, i = l,...,n, 
Ki 



(14) 
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very short time scales^ and hence is relevant for online optimization and fine-tuning under varying 
operating conditions that affect 0L2. Since equation ([12]) can be satisfied exactly, this approach leads 
to the fastest possible asymptotic convergence rate to the synchronous state among all systems with 
the same network structure and identical /3j values. 

More important, even when the /3j values are not identical, as in the case of real power grids. 



adjusting Ri according to equation ( [13] ) or Di according to equation ( ]T4] ) in general results in sub- 
stantially improved stability. For example, consider the three test systems and the Guatemala power 
grid in Table 1, which are all stable before any parameter adjustment. As these systems have non- 
identical /3j values (Table 1, column 7), their stability was determined directly from equations ([6]) 
and ([7]) before diagonalization; the stability is quantified by X^^^, the largest Lyapunov exponent 
excluding the null exponent associated with the uniform shift of phases. A naive approach based 
on parameter adjustment that homogenizes /3j to their average /3 = does not necessarily 

improve stability, as illustrated by the opposite effect observed in the Guatemala system (Table 1, 
columns 8 and 10). However, the adjustment to /3opt improves stability in all systems, by a factor 
ranging from 5 to nearly 90 in terms of }^^^^ (Table 1, column 11). 

This demonstrates that the mechanism by which synchronization stability can be enhanced 
requires not only homogenization of the function /3j of the generator parameters ifj, D-i, and 
but also matching of this function with the structural and the dynamical state of the network rep- 
resented in 0L2- Furthermore, we can show that the maximum Lyapunov exponent }^^^^ is locally 
minimum at /3j = /3opt along any given direction in the /3i-space for all systems we consider, which 
is a necessary condition for this parameter assignment to be locally optimal for synchronization 
stability. The approach is therefore appropriate for enhancing the stability of power networks for 
which the synchronous state is already stable. If the synchronous state is not stable, as in the case 
of the Northern Italy and Poland systems for the dynamic parameters we consider, we can first 
make it stable by adjusting the generators' transient reactances x'^ ^ (see Methods and Table 1, col- 
umn 9). Then, the synchronization stability can be further enhanced by the adjustment of /3j to 
/3opt, as illustrated in Fig. [3] for the Northern Italy system. Taken together, these results provide a 
systematic approach to strengthen stability in power- grid networks. 



Discussion 

Power grids are dynamic entities whose structure is history-dependent and evolves in a decentral- 
ized way that is often determined by market. The consequent difficulty to rationally design and 
modify the network may be seen as a major limiting factor in optimizing synchronization in these 
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systems. Previous theoretical studies of synchronization in oscillator networks have shown that the 
structure of the interaction network is a determinant factor for the dynamical units to synchronized. 
In the most studied case of diffusively-coupled oscillators, it has been further shown that certain 
network structures that inhibit synchronization can in fact facilitate synchronization when in pres- 
ence of other synchronization-inhibiting network structures, such as negative interaction^. But in 
power grids, the relevant network is not simply the physical network of transmission lines, and this 
causes the structure and dynamics to be more intimately related than in such idealized models (see 
Methods). In part because of this, here we have been able to show that the stability of synchronous 
states in power grids can be enhanced by tuning parameters of the dynamical units rather than the 
network. 

Such enhancement can be implemented through the condition expressed by (3i = 2^/1x2. The 
r.h.s. of this equation accounts for the network structure and indirectly depends on the generators 
(see Methods) while the l.h.s. depends only on the generator parameters Hi, Ri, and Di, among 
which the latter two do not appear on the r.h.s. Of these, the droop parameter Ri offers a solution 
that accounts for slow changes in demand and the long-term evolution of the network. We argue 
that the other parameter, the damping coefficient Di, has the flexibility required to cope with rapid 
changes caused by a fault or fluctuation, since the adjustment of this parameter can be realized 
through very fast control loops (e.g., by adding power system stabilizers). Our optimization scheme 
can be effective even when modeling load dynamicd^ (Fig. |4]), and it is complementary to an 
approach proposed recentlj^^to mitigate saddle-node instability by adjusting power scheduling or 
line impedances (Fig. [5]). Conversely, the adjustment of power scheduling and/or line impedances 
can also be effective within the reduced network model used in our study (Fig. |6]). 

We suggest that our findings are potentially important for ongoing research on smart grids, 
which is making it ever more important to understand optimization and control of power-grid 
dynamics. While the study of stabilization via damping of oscillations has a long historjl^^I^, 
the proposed coordinated tuning of generators to enhance stability is a timely approach in view 
of the upcoming availability of system- wide data from phasor measurement unitP^. These data 
will provide accurate information about the synchronization state of the power grid, which can 
be integrated with the online control of generator parameters. Therefore, besides contributing to 
the devise of more robust systems, our findings provide insights for the development of efficient 
controllers. Such controllers may help advance research on self-healing systems that can recover 
from failures in real time. This can lead to systems that are less prone to cascading failures, 
which, as amply discussed in the literature on interdependent network^, have consequences that 
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transcend the power grid itself. 

Our analysis of power-grid dynamics also provides a fresh view of synchronization in com- 
plex networks in general. By identifying as the leading factor for stability the relation between 
the specifics of the dynamical units and the network structure, it contrasts with most previous 
studies, which focused on the role of the network structure alone (analyses of correlations be- 
tween the natural frequencies of phase oscillators and their connection topology are among the 
few exceptiond^^^. Aside from power grids and other technological applications, this is likely 
to have implications for natural systems, such as biological ones, where the dynamical units co- 
evolve with the network structure!^. Moreover, our analysis establishes a master stability condition 
for the synchronization of power generators, casting the problem in the framework of oscillator 
networkd^ used to investigate collective behavior and nonlinear phenomena in a wide range of 
complex systems. This facilitates comparison across different domains, and, we hope, will inspire 
similar developments on the tailoring of collective behavior in other real systems. 

Methods 

Power flow calculations and synchronous states. The power and voltage at each node were 
determined given the net injected real power and the voltage magnitude at each generator node, 
the power demand at each non-generator node, the reference voltage phase, the admittance of each 
power line, and the capacities of the network components. These power flow calculations were 
performed using the Power Systems Analysis Toolbox (PSATj^. To identify the synchronous 
states, the phase 6* of generator i was determined from the power and voltage at node i assuming 
the classical modeP°. In this model, the generator is represented by a voltage source with constant 
magnitude Ei and variable phase 5i that is connected to the rest of the network through a transient 
reactance x'^ ,■ and a terminal. Throughout the paper, we represent each generator and its terminal 
by a single node, except in the Kron reduction, where the terminals are treated as independent 
nodes to be eliminated. 



Derivation of the swing equation. Equation ([2]) can be derived by setting the rate of change of 
the angular momentum of the rotor equal to the net torque acting on the rotor: 

where J is the moment of inertia in kg-m^, T^i is the mechanical torque in N-m accelerating the 
rotor, and Tgj is the typically decelerating torque in N-m due to electrical load in the network. 
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Multiplying both sides by and using the fact that the torque in N-m multiplied by the angular 
velocity in radians per second gives the power in watts, the equation can be written in terms of 
power: 

cPSi 

JuJi— = Pmi- Pet- (16) 

To make Pmi and Pet per unit quantities, we divide both sides of the equation by the rated power Pr 
(used as a reference). The factor Juji then becomes 2Hi/ijji, where we defined the inertia constant 
Hi = Wi/Pji (in seconds) and the kinetic energy of the rotor Wi = \Jujj (in joules). Noting that 
oji is approximately equal to the reference frequency ujr in systems close to synchronization, we 
obtain equation ([2]). A more detailed description of this derivation can be found in ref . [30] (second 
edition, pp. 13-16). 



Relation to the master stability formalism. Equations (|6]) and Q when written individually for 
each node i, can be expressed as 
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This is in the same general form as the variational equation for the class of coupled oscillators 
considered in ref. 34 which for a synchronous state s = s(t) is 

d 



dt 



(18) 



where Xj = F(xi) describes the node dynamics, H is the coupling function, and aGij represents 



the strength of coupling from node j to i. The factors DF{s) and DH(s) in equation ( [18] ) are both 



constant matrices in equation ( [17] ), and Pij in equation ( [17] ) corresponds to crGij in equation ( [18] ), 
which are relations that can be used to derive A^(a) from equations (|6]) and ([7]) based on the results 
in ref. 34 Therefore, while the formalism in ref. 34 cannot be directly applied to equation ([2]), it 
can be applied near the synchronization manifold, which is a procedure previously proposed for a 
broad class of coupled oscillator^. 



Adjusting transient reactances for stability. We first observe that for Bij > 0, which is the 
most frequent case (Table 4), the cosine term in equation ([8]) is a destabilizing factor if \5*j\ > 
n/2. This instability often arises when a few specific generators have phases that are significantly 
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different from those of the other generators in the synchronous state, which can be due to these 
generator's transient reactance x'^^. By identifying these generators through the negative sign of 
the corresponding diagonal elements of P and reducing their x'^ ■, the instability can be suppressed. 
Applying this procedure to the Northern Italy and Poland networks indeed turns the systems stable 
(Table 1, column 9). 



Relation between structure and dynamics. The interactions between generators in a power grid 
are in principle determined by the effective admittances, which are dominantly imaginary (Table 2). 
Thus, these interactions too have a sign; an inductive reactance (Bij > 0) suggests a positive 
interaction, and a capacitive reactance (Bij < 0) a negative one. However, it can be seen from 
equation ([8]) that the interactions between the generators are in fact determined by Bij cos6*j. In 
order for 0^2 > 0, as required for stable synchronization, the terms B^ cos 5*j have to be dominantly 
positive, as observed in real systems (Table 4). One can make any of these terms positive by 
having B^j and cos 5*j both positive or both negative, where having both factors positive is far 
more common in real systems than both negative (Table 4). While Bij is mainly structural given 
the power flow solution, the parameter 5*j can be adjusted by changing specific properties of the 
generators, such as their transient reactances. Therefore, the interactions are in reality determined 
not only by the network of effective admittances, but also by the system-level (alternate current) 
dynamics, which is in turn molded by the properties of the dynamical units. As a result, important 
aspects in the network structure can be emulated by modifying tunable parameters of the dynamical 
units. 



Power-grid data. The data required for power flow calculations were obtained as follows. For 
the 10-generator system, known as the New England test system, the parameters were taken from 



ref. B91 For the 3- and 50-generator systems, the parameters were taken from refs. 30 and 50 



respectively. The data for the Guatemala and Northern Italy systems were provided by F. Milano 
(University of Castilla - La Mancha), and the data for the Poland system were provided as part 
of the MATPOWER software package!^. The dynamic data required to calculate the synchronous 
state, determine its stability, and simulate equation (|2]), are not all available for the real power 
grids, and were obtained as follows. For all systems, we assumed that before any optimization 
the damping coefficient and droop parameter satisfy Di + 1/Ri = 50 per unit for all generators. 
The parameters x'^ ^ and Hi are available for the three test systems from the respective references 
mentioned above. For the Guatemala, Northern Italy, and Poland systems, we estimated x'^,- and 
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Hi using the strong correlation observed in the test systems between each of these parameters and 
the power Pj injected by generator i into the network, as shown in Fig. |7} The estimated values 
are x'^^ ~ 92.8-F^~^ '^ and Hi ^ 0.04Pj, where Pi is in megawatts, x'^i is in per unit, and Hi is in 
seconds. 
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Figure 1: Physical versus effective network for the power grid of Northern Italy, a, Rep- 
resentation of the physical network of transmission lines, which has 678 nodes and 822 links, 
b, Representation of the network of effective admittances, which is an all- to-all network with 170 
nodes corresponding to the generators in the system (a subset of the nodes in panel a). The color 
scale of the lines indicates the link weights, ranging from yellow to red to black (scaled differ- 
ently for each panel), defined as the absolute value of the corresponding admittance. For clarity, in 
panel b we show only the top 50% highest- weight links. 
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Figure 2: Stability of synchronous states for systems with (3i = (3. a, Real part of the Lya- 
punov exponents A_|_ (continuous lines) and A_ (dashed lines) as functions of a for increasing 
values of (3 in the range 1.0-3.0. Although changing (3 causes changes in the shapes of ReX±, the 
region of stability defined by A/3(a) < is always a > 0. b, Real part of the Lyapunov exponent 
A+ as a function of (3, which has a minimum at /3 = 2y/a (illustrated here for a = 1). The stability 
of synchronous states is maximum at this point. 
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Figure 3: Enhancement of the stability of synchronous states, a-d, Response of the original 
network (a), the network with adjusted x'^ ■ (b), the network with A = /3 (c), and the network with 
A = /5opt (d) to a perturbation in the Northern Italy power grid. The perturbation was applied to 
the phase of each generator in the synchronous state at t = 0, and was drawn from the Gaussian 
distribution with mean zero and standard deviation 0.01 rad. In each case, the network layout is the 
same as in Fig. lb and the bottom panel shows the time evolution of ^-^ = 5ii — 5*^, where Sn is the 
phase of generator i relative to generator 1 and 5*^ is the corresponding phase in the synchronous 
state. Generator 1, shown as a white node in the network, is used as a reference to discount phase 
drifts common to all generators. The other nodes and their time-evolution curves are color-coded 
by the maximum value of \S^i \ for 2 < t < 3. By adjusting the transient reactance of the generators, 
the divergence from the unstable steady state is converted to exponential convergence (a and b). 
This stability is improved upon adjusting the generator parameters to ensure a common value for 
(3i, but is further improved when this common value is tuned to /3opt = 2^/02 (c and d). 
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Figure 4: Enhancement of synchronization stability in power-grid model incorporating load 
dynamics. Colors indicate the factor by which the stability improves after adjusting all (3i to the 
value /3opt, where the darkest red corresponds to any factor > 7 and black to any factor < 1. 
The stability of a synchronous state is measured by the largest nonzero Lyapunov exponent X^^^ 
computed within the Bergen-Hill model'^'. In this model, instead of eliminating non-generator 
nodes by the Kron reduction, load dynamics is modeled by assuming that the real power is a linear 
function of the voltage frequency for all power-consuming nodes, while the voltage magnitude is 
assumed to be constant for all nodes. For illustration, we use the 3-generator system studied in ref. 



48 with the power injection modified to Pi = 5,6,7 per unit. The improvement factor is shown as 



a function of the parameters Dg^n and -Dioad, where Dg^n is a coefficient given by {Di + 1/ Ri)/ujR 
(assumed be the same for all generator nodes) and -Dioad is the frequency coefficient (assumed 
to be the same for all the other nodes). In most cases, our method yields significant stability 
improvement, demonstrating its effectiveness beyond the reduced network model. 
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Figure 5: Combination of complementary approaches to mitigate instabilities associated 
with saddle-node bifurcations. As an example, we use the same 3-generator system described in 
the caption of Fig. |4] with power injection Pi = 7.994, 3.006, 7.000, simulated using the Bergen- 
Hill model^. For Dge^ = -Dioad = 1 per unit, the system is near a saddle-node bifurcation, and 
thus the Lyapunov exponent X^^"^ is negative real and close to zero. At each iteration of the gradi- 
ent descent-like method described in ref. 48^ which changes the power injected by the generators 
(independently of the values of /3j and thus of Di), we compute A™^"^ both with and without ad- 
justing Pi to /3opt by tuning the damping coefficients Di. The improvement factor, measured by 
the ratio between the smallest \^^^ obtained through the iterative process with and without the 
adjustment, is shown as a function of -Dgen for different values of -Dioad- In most cases, significant 
additional improvement results from the adjustment of (3i, illustrating that near instabilities the two 
methods can be combined to achieve enhancement not possible by either method alone. This re- 
sult is robust against heterogeneity in the network parameters, as illustrated in the inset histogram 
for the system with Dgen = 0.5 and -Dioad = 1 (blue dot in the main plot), where each of these 
coefficients is independently perturbed according to the Gaussian distribution with mean zero and 
standard deviation 0.1. 
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Figure 6: Enhancement of synchronization stability via tuning of power injection within 
the reduced network model. We use the same 3-generator system as in Figs. |4] and |5} where 



all system parameters are the same as in ref. 48 except for the power injection and consumption, 



which are initially set to 60% of the values used in that reference. This makes the system unstable, 
as indicated by the positivity of the Lyapunov exponent A™*^"^ for the synchronous state. The power 
injections were then adjusted iteratively according to the gradient descent-like scheme of ref. |48| 



and X^^"^ is plotted against the number of iterations for several values of the coefficient -Dgen. In all 
cases, the power adjustment stabilizes the synchronous state, which is a consequence of increasing 
a2 across zero, as shown in the inset. Note that a2 is independent of Dgen because the synchronous 
state in the reduced network model is unaffected by changes in -Dgen- 
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Figure 7: Relationship between the power produced by individual generators and their pa- 
rameters in the three test systems of Table 1. a, The transient reactance x'^ ^ versus the power 
Pi injected into the network, b, The inertia constant versus Pj. We used the approximate func- 
tional relations revealed by these data to estimate the transient reactance and the inertia constant 
for each generator in the Guatemala, Northern Italy, and Poland systems. 
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Table 1 : Structural properties and synchronization stability of the systems considered. 



System heterogeneityt Synchronization stability* 



System§ 


Nodes 


Linl^s 


Generators 


Yo 


Y 


ft 


Original 


x'^ ^ adjusted 


ft = ,s 


Pi = 0opi 


3-generator test system 


9 


9 


3 


0.39 


0.09 


0.83 


-1.71 




-2.21 


-8.69 


10-generator test system 


39 


46 


10 


0.81 


0.38 


0.37 


-0.24 




-0.37 


-3.65 


50-generator test system 


145 


422 


50 


1.91 


1.21 


2.07 


-0.02 




-1.53 


-1.75 


Guatemala power grid 


370 


392 


94 


6.86 


1.18 


1.22 


-0.33 




-0.02 


-4.19 


Northern Italy power grid 


678 


822 


170 


3.42 


0.84 


2.02 


7.20 


-0.48 


-1.59 


-3.70 


Poland power grid 


2383 


2886 


327 


2.26 


3.00 


0.92 


140.07 


-0.03 


-0.02 


-1.53 



tAs a measure of heterogeneity, we show the standard deviation normalized by the average. The quantities considered are the weighted degrees 
computed for Yo and Y, which represent the structure of the physical and effective networl<, respectively, and the parameter ft, which represents 
properties of the generators. ^As a measure of synchronization stability, we show the Lyapunov exponent A™*™. We consider this exponent for the 
original parameter values, for the transient reactance x'^ ^ adjusted, for all ft set equal to their average and for the generator parameter and/or 
Di adjusted to ensure ft = ftpt. §See Methods for a description of data sources. 
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Table 2: Phase difference in synclironous states and tlie real versus imaginary 
parts of the admittances in the physical and effective networks. 



Physical network^ Effective network 



System 


Mean |6*.|/7r 


Mean \Gnn\ 


Mean IBnul 


Mean \Gij\ 


Mean |-Bij| 


3-generator system 


0.06 


0.95 


12.0 


0.2367 


1.275 


10-generator system 


0.07 


6.81 


80.9 


0.3784 


1.036 


50-generator system 


0.12 


2.57 


46.3 


0.1978 


1.242 


Guatemala 


0.08 


7.74 


62.1 


0.0015 


0.004 


Northern Italy: original 


0.07 


105.88 


758.3 


0.0106 


0.039 


xjj ^ adjusted 


0.05 






0.0104 


0.042 


Poland: original 


0.25 


24.06 


597.3 


0.0012 


0.018 


x'^ ^ adjusted 


0.13 






0.0013 


0.019 



tThe real and imaginary components of the admittance -ioij are denoted -Goij and -Bo<j . respectively. 



Table 3: 2-norm of the symmetric and antisym- 
metric parts of the matrix P''. 



System 


Symmetric 


Antisymmetric 


3-generator system 


0.68 


0.011 


1 0-generator system 


0.48 


0.018 


50-generator system 


1.89 


0.016 


Guatemala 


3.04 


0.085 


Northern Italy; original 


2.45 


0.030 


x'^ ^ adjusted 


34.27 


0.068 


Poland: original 


14.45 


0.261 


x'j^ ^ adjusted 


234.24 


0.982 
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Table 4: Fraction of positive off-diagonai eiements in 
tlie matrices (Bij), (cos5* ), and (Bij cosS*-). 



System 




cos (5,*,. 


Bij cos Sf, 


3-gGnGrator systGm 


1.000 


1 .000 


1 .000 


1 0-generator system 


1.000 


1.000 


1.000 


50-generator system 


0.689 


0.995 


0.687 


Guatemala 


0.999 


1.000 


0.999 


Northern Italy; original 


0.999 


0.973 


0.973 


x'^ ^ adjusted 


0.997 


1.000 


0.997 


Poland: original 


0.978 


0.829 


0.810 


x'j^ ^ adjusted 


0.957 


0.992 


0.949 
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